#include "simulatorParameters.h"

SimulatorParameters::SimulatorParameters()
{
	TimeStepCounter = 0;
}

SimulatorParameters::~SimulatorParameters()
{
	MIter_RockProperties mit = mapRockProp.begin();
	for (; mit != mapRockProp.end(); mit++)
	{
		double *K = mit->second->K;
		delete[] K; K = 0;
	}
	mapRockProp.clear();
	mapBC.clear();	
}
	
void SimulatorParameters::loadParameters()
{
	int flag;
	double val;
	string str[6];

	ifstream fid;
	fid.open("datafiles/simulator-parameters");
	if ( !fid.is_open() )	
		throw Exception(__LINE__,__FILE__,"flag-list not found or it does not exist\n");

	// start reading file
	// -------------------------------------------------------------------------
	setPositionToRead(fid,"# C.F.L - AppOrder - numTimeSteps - totalSimulationTime - EF - RF:");	

	fid >> str[0] >> str[1] >> str[2] >> str[3] >> str[4];
	_CFL = strtod(str[0].c_str(), 0);
	_AppOrder = atoi(str[1].c_str());
	_NOTS = (!str[2].compare("*"))?0:atoi(str[2].c_str());
	_TTOS = (!str[3].compare("*"))?0:atoi(str[3].c_str());		
	_EF = atoi(str[4].c_str());
	_RF = atoi(str[5].c_str());
	
	setPositionToRead(fid,"# abstol - rtol - stol - maxit - maxf:");
	fid >> str[0] >> str[1] >> str[2] >> str[3] >> str[4];
	_abstol = strtod(str[0].c_str(), 0);
	_rtol = strtod(str[1].c_str(), 0);
	_stol = strtod(str[2].c_str(), 0);
	_maxit = atoi(str[3].c_str());
	_maxf = atoi(str[4].c_str());
	
	setPositionToRead(fid,"# Oil density, Water density, Oil viscosity"
			", Water viscosity:");
	fid >> str[0] >> str[1] >> str[2] >> str[3];
	_oil_density = strtod(str[0].c_str(), 0);
	_water_density = strtod(str[1].c_str(), 0);
	_oil_viscosity = strtod(str[2].c_str(), 0);
	_water_viscosity = strtod(str[3].c_str(), 0);	
	
	setPositionToRead(fid,"# Mesh dimension:");
	fid >> str[0];
	_meshdim = atoi(str[0].c_str());
	
	setPositionToRead(fid,"# dom, K and phi:");
	readRockProperties(fid);
	
	setPositionToRead(fid,"# Dirichlet:");
	mapFlag(fid,"dirichlet");
	
	setPositionToRead(fid,"# Neumann:");
	mapFlag(fid,"neumann");
	
	setPositionToRead(fid,"# Saturation:");
	mapFlag(fid,"saturation");
	
	setPositionToRead(fid,"# Initial conditions: water saturation:");
	fid >> str[0];
	_Sw_initial = strtod(str[0].c_str(), 0);
}

// skip comments or any other unnecessary string until a specific point be reached 
// =========================================================================================
// IN:
// fid - pointer to read data from file
// str2 - file is read until str2 is found
// OUT:
// fid - ready to read real data
void SimulatorParameters::setPositionToRead(ifstream &fid, string str2)
{
	string str1;
	do{
		getline(fid,str1);
	}while( str1.compare(str2) );	
}

void SimulatorParameters::readRockProperties(ifstream &fid)
{	
	int i;
	const int dim = meshDimension();
	RockProperties *rockprop;
	string str;
	
	while( true )
	{
		fid >> str;
		if (!str.compare("end")) break;
		
		rockprop = new RockProperties;
		
		// get domain flag
		const int dom = atoi(str.c_str());
		
		// get permeability tensor
		rockprop->K = new double[dim*dim];
		for (int i=0; i<dim*dim; i++)
		{
			fid >> str;
			rockprop->K[i] = strtod(str.c_str(), 0);
		}
		
		// get porosity
		fid >> str;
		rockprop->porosity = strtod(str.c_str(), 0);
		
		mapRockProp[dom] = rockprop;
	}	
	if (mapRockProp.size()==0)
		throw Exception(__LINE__,__FILE__,"You must provide rock properties\n");
	else
		printf("number of dimension provided: %d\n",mapRockProp.size());
}

// store flags and their values into a specific map container
// =========================================================================================
// IN:
// fid - pointer to read data from file
// OUT:
// mapflag - a map container: flag -> value
void SimulatorParameters::mapFlag(ifstream &fid, string whatmap)
{
	string str;
	BdryConditionType *bct;

	while( true )
	{
		fid >> str;
		if (!str.compare("end")) break;
		
		int flag = atoi(str.c_str());

		fid >> str;
		double val = strtod(str.c_str(), 0);

		// for every flag number create a new BdryConditionType pointer and add
		// inform type of condition and its value.
		bct = new BdryConditionType;
		bct->type = whatmap;
		bct->val = val;

		// two map container are used: one for saturatio and other to dirichlet/neumann
		// boundary conditions
		if (!whatmap.compare("saturation"))
			mapSaturation[flag] = bct;
		else
			mapBC[flag] = bct;
	}	
}

const double* SimulatorParameters::getPermeability(const int &dom) 
{
	MIter_RockProperties mit = mapRockProp.find(dom);	
	if (mit != mapRockProp.end())
		return mit->second->K;	
		
	cout << "Warning: getPermeability() return null value\n";
	cout << __FILE__ << "\t at line " << __LINE__ << endl;
	return 0;
}
	
double SimulatorParameters::getBC_Value(const int &flag)
{	
	MapFlagIter mIter = mapBC.find(flag);
	if (mIter != mapBC.end())
		return mIter->second->val;
	
	cout << "Warning: getBC_Value() return null value\n";
	cout << __FILE__ << "\t at line " << __LINE__ << endl;
	return 0;	
}

bool SimulatorParameters::hasNodeWell(const int &flag)
{
	// if node has a prescribed saturation it's supposed that exist a well	
	if ( mapSaturation.size() )
	{
		MapFlagIter MIter = mapSaturation.find(flag);
		return (MIter != mapSaturation.end());
	}
	else
	{
		isNodeFree(flag);
		return well;	
	}
}

bool SimulatorParameters::isNodeFree(const int &flag)
{
	MapFlagIter mIter = mapBC.find(flag);
	if (mIter == mapBC.end())
		return true;
	else
	{
		BdryConditionType *bct = mIter->second;
		return ( !bct->type.compare("dirichlet") )?false:true;
	}
}

bool SimulatorParameters::finishSimulation(double t)
{
	if ( numTimeStep() )
		return (++TimeStepCounter > numTimeStep())?true:false;

	return (t > (double)totalTimeSimulation())?true:false;	
}

// chose between first order and second order approximation 
bool SimulatorParameters::isHigherOrder() const
{
	return _AppOrder;
}

// numeic parameters
double SimulatorParameters::CFL() const
{
	return _CFL;
}

int SimulatorParameters::numTimeStep() const
{
	return _NOTS;
}

int SimulatorParameters::totalTimeSimulation() const
{
	return _TTOS;
}

int SimulatorParameters::EF() const
{
	return _EF;
}

int SimulatorParameters::RF() const
{
	return _RF;
}	

double SimulatorParameters::abstol() const
{
	return _abstol;
}

double SimulatorParameters::rtol() const
{
	return _rtol;
}

double SimulatorParameters::stol() const
{
	return _stol;
}

int SimulatorParameters::maxit() const
{
	return _maxit;
}

int SimulatorParameters::maxf() const
{
	return _maxf;
}	

// physical parameters
double SimulatorParameters::getPorosity(const int &dom)
{
	MIter_RockProperties mit = mapRockProp.find(dom);	
	if (mit != mapRockProp.end())
		return mit->second->porosity;
		
	cout << "Warning: porosity() return null value.";
	cout << "Domain "<< dom << " not found.\n";
	cout << __FILE__ << "\t at line " << __LINE__ << endl;
	
	return 0;
}

double SimulatorParameters::oilDensity() const
{
	return _oil_density;
}

double SimulatorParameters::waterDensity() const
{
	return _water_density;
}

double SimulatorParameters::oilViscosity() const
{
	return _oil_viscosity;
}

double SimulatorParameters::waterViscosity() const
{
	return _water_viscosity;
}	

double SimulatorParameters::Sw_initial() const
{
	return _Sw_initial;
}	

int SimulatorParameters::meshDimension() const
{
	return _meshdim;
}

// DEBUG PURPORSES
void SimulatorParameters::printParameters() //printFlags
{
	printf("%f %d %d %d %d %d\n",
			CFL(),_AppOrder,numTimeStep(),totalTimeSimulation(),EF(),RF());	
	printf("%f %f %f %d %d\n",
			abstol(),rtol(),stol(),maxit(),maxf());
	printf("%f %f %f %f\n",
				oilViscosity(),waterViscosity(),oilDensity(),waterDensity());
	printf("%d \n",meshDimension());
	
	int dim = meshDimension();
	for (MIter_RockProperties mit = mapRockProp.begin(); mit != mapRockProp.end(); mit++)
	{
		RockProperties *rprop = mit->second;
		printf("dom[%d]\n",mit->first);
		printf("porosity: %f\n",rprop->porosity);
		printf("permeability: \n");
		for (int i=0; i<dim*dim; i++) printf("%.4E ",rprop->K[i]);
		printf("\n");
	}
	
	for (MapFlagIter mit = mapBC.begin(); mit != mapBC.end(); mit++)
	{
		BdryConditionType *bc = mit->second;
		printf("type of BC: %s\t",bc->type.c_str());
		printf("flag: %d\t",mit->first);
		printf("value: %f\t",bc->val);
		printf("\n");
	}	
	
	printf("Sw_initial: %f\n",Sw_initial());
	
}
